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Abstract. We use continuous-time quantum Monte Carlo simulations to study 
retardation effects in the metallic, quarter-filled Holstein model in one dimension. 
^ ' Based on results which include the one- and two-particle spectral functions as well as 

^ . the optical conductivity, we conclude that with increasing phonon frequency the ground 

\ state evolves from one with dominant diagonal order — 2k-p charge correlations — to one 

with dominant off-diagonal fluctuations, namely s-wavc pairing correlations. In the 
parameter range of this crossover, our numerical results support the existence of a 
spin gap for all phonon frequencies. The crossover can hence be interpreted in terms 
of preformed pairs corresponding to bipolarons, which are essentially localised in the 
Peierls phase, and "condense" with increasing phonon frequency to generate dominant 
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' pairing correlations. 

> 

^ ! 1. Introduction 

00 

The concept of preformed fermion pairs or bosonic degrees of freedom which condense 
to form a superfluid can be found in many domains of correlated quantum many-body 
C\| ■ systems. Examples include the resonating valence bond theory of high-temperature 
superconductivity PQ, Mott metal- insulator transitions in cold atoms [2], or transitions 
between charge- density- wave and superconducting states in the family of dichalcogenides 
[3]. In this article, we show in the framework of the Holstein model that these concepts 
can be carried over to one dimension where, in the absence of continuous symmetry 
breaking, the phase transition is replaced by a crossover. 

Quantum lattice fluctuations play a crucial role in one dimension. For classical 
phonons (i.e., in the adiabatic limit), the Peierls instability leads to a 2kp periodic lattice 
deformation, resulting in corresponding charge ordering of electron pairs (bipolarons) 
and the opening of a gap at the Fermi level. In contrast, for the Holstein model, 
one can argue and give numerical evidence [4] that quantum phonons will close the 
charge gap but preserve the spin gap, thereby giving rise to a "bad metal" state. In 
the latter, electrons bind into pairs and are mostly localised in the 2kp charge order 
pattern. This results in a very high effective charge-carrier mass or, equivalently, in 
a very small Drude weight. We refer to this state as "the Peierls state", which is to 
be distinguished from an insulating Peierls state with long-range charge order [5]. The 
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latter has been studied intensely in Holstein-type models with and without spin, see 
discussion in section [2j The metallic Peierls state may be regarded as a more general 
concept for quasi-one- dimensional systems, applicable to noncommensurate fillings or 
the normal state (T > T c ) of a Peierls insulator [6]. On the level of bosonization, 
a metallic state with dominant 2/cp charge correlations requires a finite spin gap [6]. 
In the anti-adiabatic limit, the lattice adjusts instantaneously to the electronic motion, 
and the Holstein model maps onto the attractive Hubbard model with dominant pairing 
correlations. 

The question addressed in this article is the nature of the crossover from the 
adiabatic to the anti-adiabatic regime. The limiting cases of low (but finite) and 
high (or infinite) phonon frequency can be described as a Luther-Emery liquid [7J [8] , 
corresponding to a Luttinger liquid with a gapless charge mode and a gapped spin mode. 
A key result of the present work is that the crossover occurs without a closing of the spin 
gap, which suggests that the bipolarons existing in the Peierls state do not dissociate 
as a function of phonon frequency but instead constitute the preformed pairs of the 
superfluid phase. 

The article is organised as follows. In section |2] we introduce the model and explain 
the key ideas of the continuous-time quantum Monte Carlo method. Numerical results 
are presented in section [31 and we end with conclusions and a discussion of open issues 
and future research directions in section HI 



2. Model and method 



The one- dimensional Holstein model is defined by the Hamiltonian [9] 
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Qi (n% - i) , 



with the tight-binding dispersion relation = — It cos(ka)—fi and the chemical potential 
fi. The operator c| CT creates an electron in the Wannier state centred on lattice site i 
with spin a, c\ a = L~ 1 / 2 ^ . e lfcj, cj CT creates an electron in a Bloch state with wavevector 
k and spin a, and hi = c\ a c ia is the particle number operator; Qi and Pi denote the 
displacement and momentum of the harmonic oscillator at site i. 

To solve this model numerically without approximations (except for the finite 
system size), we integrate out the phonons at the expense of a retarded density-density 
interaction. The partition function of the resulting model is given by the path integral 

[dcW] e - {So+s ^ ] (2) 

with the following contributions to the action: 

So = f dr 4a(r) " ^) c ja (r) , (3) 
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where ijj = 1 for nearest neighbours and zero else, and 

S cp = f dr [ dr' - 1] D {i - j,r - r') [n {r') - 1] . (4) 

Jo Jo id 

Here, c\ a {r) is a Grassmann variable, (5 denotes the inverse temperature, and D (i — 
j, t — t') is the phonon propagator. For dispersionless Einstein modes, D Q is diagonal 
in space and takes the form 

D (z-j,t-t') = 5^P(t-t>), (5) 

= uoe +e 

v i 2 1 - e-#"» 

The phonon frequency is given by ujq = a/ K/M . In the following, we use t as the unit 
of energy, and set h and the lattice constant to unity 

In the fermionic form given by equations ([3]) and fll]), the problem can be solved 
very efficiently by means of interaction-expansion continuous-time quantum Monte Carlo 
methods pUl ELlJ E2] • The advantage of such a formulation is that we do not have to 
sample the lattice degrees of freedom directly, thereby eliminating the need for a Hilbert 
space cutoff. The details of the implementation of this algorithm can be found in [13] . 
The CTQMC method has previously been applied to the spinless [H] and the spinful 
Holstein model [15], as well as to a model with nonlocal electron-phonon interaction 
|16j . Simulations are carried out on one-dimensional lattices of length L using periodic 
boundary conditions. 

The one-dimensional, half-filled Holstein model — with or without an additional 
Hubbard term to describe electron-electron repulsion — has been studied for many 
years. The existence of a metal-insulator transition at zero temperature has been 
established with the help of exact numerical methods, including quantum Monte Carlo 
PH UHl HH 1201 US] the density-matrix renormalization group @J M\ 1221 [231 123 [25] and 
exact diagonalization [26]. An overview of analytical work on the Peierls transition can 
be found in [27] . Many aspects of the relevant physics can also be captured by a spinless 
model with filling n = 0.5, including an extended metallic region in the phase diagram 
[28| |29| |30| |3~T] . the renormalization of the phonon mode [301 E21 ESI El], and soliton 
excitations [14]. Retardation effects in the half-filled Holstein model were discussed, for 
example, in [351 ESj, and for the spin-Peierls transition in [37] . 

Here we study the effect of phonon frequency in the quarter-filled Holstein model 
at a fixed electron-phonon coupling strength. For quarter filling and low phonon 
frequencies, the absence of first-order umklapp scattering — dominant at half filling — 
suppresses the Peierls instability, and stabilises a state with gapless charge but gapped 
spin degrees of freedom, as well as a nonzero single-particle excitation gap, and dominant 
but critical 2kp charge correlations (a Luther-Emery liquid with interaction parameter 
K p < 1). The spin gap is a result of backscattering. Signatures of this Luther- 
Emery phase have also been observed in the Holstein-Hubbard model with comparable 
electron-phonon and electron-electron interactions [381 LTHl [191 E21 [23]. Our results 
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presented below suggest that in contrast to half filling, we can observe dominant pairing 
correlations for high phonon frequencies at a commensurate density; dominant pairing 
for incommensurate filling has been reported before [39]. Recently, the existence of a 
correlated singlet state in the Holstein-Hubbard model away from half filling has been 
suggested jlQl ST]. Compared to half filling, the critical value for the transition to 
a Peierls insulator is significantly larger at quarter filling [19]. This increase may be 
understood in terms of a devil's staircase, similar to charge- density- wave transitions in 
extended Hubbard models away from half filling [8]. As for half filling, the critical value 
of the electron-phonon coupling strength for the metal-insulator transition increases 
with increasing phonon frequency as a result of enhanced lattice fluctuations |39j . 



3. Results 



In the anti-adiabatic limit ujq — » oo, P(r) in equation (jSJ) reduces to a Dirac 5-function, 
facilitating the above-mentioned mapping of the Holstein model onto the attractive 
Hubbard model, with U = g 2 /k. The ratio of this binding energy and the bandwidth 
W — At gives the dimensionless electron-phonon coupling 



Throughout this article, we set A = 0.35 and concentrate on the quarter-filled band with 
k-p = 7r/4. In the following, we first establish a picture of the physics on the basis of 
equal-time correlation functions, and then turn to dynamical correlations. 



3.1. Static correlation functions 

Figure [TJ shows equal-time correlation functions for charge, spin, and pairing, as well as 
the momentum distribution function, at various phonon frequencies. In the adiabatic 
limit uo/t = 0, any A > leads to an insulating state. As discussed in section [21 
we choose the coupling strength A = 0.35 such that we have a metallic Luther-Emery 
liquid with dominant 2/cp charge correlations for low phonon frequencies, and then study 
the evolution as a function of increasing uo/t. In particular, we have verified that for 
LJo/t = 0.1, the lowest nonzero frequency considered in the following, there is no long- 
range order; this can be seen from the finite-size dependence of the charge susceptibility 

The density (or charge) structure factor, defined as 

iV(g) = ^e^((n r n )-(n r )(n )), (7) 

r 

is plotted in figure QJa). For classical phonons, u = 0, the Peierls instability leads to 
long-range 2/cf charge order at zero temperature. As discussed above, quantum lattice 
fluctuations (occurring for ojq > 0) can melt this order, and lead to a state with dominant 
but power-law charge correlations [15], as confirmed by the cusp at 2kp = ir/2 in 
figure [U(a). The magnitude of the peak at q = 2kp initially decreases and then saturates 
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Figure 1. Static correlation functions for different values of the phonon frequency 
LOo/t at A = 0.35, and for a quarter-filled band (n = 0.5). The panels show (a) the 
charge structure factor, (b) the pairing correlator, (c) the spin structure factor, and 
(d) the momentum distribution function. Here L = 28 and fit = 40. 



upon increasing the phonon frequency, signalling competing ordering mechanisms as 
well as enhanced lattice fluctuations. The linear form of the charge structure factor 
at long wavelengths [see figure [T]( a)] indicates a 1/r 2 power- law decay of the real-space 
charge correlations and hence a metallic state. 

In figure [T^b), we present the pair correlation function in the onsite s-wave channel, 

P(r) = (At A >, At = 4 t c r V (8) 

In contrast to the density correlator which picks up diagonal order, P(r) detects off- 
diagonal order characteristic of a superconducting state. In the Peierls state obtained 
for classical phonons, diagonal long-range charge order leads to an exponential decay of 
pairing correlations at long distances. The fluctuations resulting from a finite phonon 
frequency close the charge gap, and render the pairing correlations critical. Comparing 
figures [T](a) and [H(b), we see that the suppression of the 2/cp charge correlations is 
accompanied by an increase of the pairing correlations, especially at large distances. 
A possible interpretation is that with increasing phonon frequency, the trapping of 
bipolarons in the 2kp lattice modulation gives way to a " condensation" (in the usual 
sense of superfluidity in one dimension) of those preformed pairs. 

The above interpretation relies on the electron pairs remaining bound upon 
increasing the phonon frequency. Evidence for the existence of bound pairs is provided 
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Figure 2. Pairing (□) and charge-density-wave (O ) correlations at the largest 
distance r = L/2 for different system sizes L. Here A = 0.35, fit — L and (a) u>o/t= 0.1, 
(b) u!o/t — 4. Lines are fits to the form f(r) = A/r 11 . 



by the equal-time spin-spin correlation function 

S(q) = ^^(S' r S' ), (9) 

r 

which we show in figure (He). For all values of the phonon frequency, the cusps at q = 
and q = 2kp are smeared out, thus lending support to an exponential decay of the spin 
correlations. The presence of a spin gap will be explicitly confirmed by results for the 
dynamical spin structure factor below. 

Finally, we plot in figure QJd) the momentum distribution function 

^) = E<4a.>- (io) 

Following the above interpretation, we expect a smooth variation of this quantity as a 
function of wavevector k. In particular we do not expect any nonanalytical behaviour at 
k = kp or k = 3kp. For classical phonons, the smooth variation of n(k) stems from the 
nonzero single-particle gap of the Peierls state. In the anti-adiabatic limit, where the 
model maps onto the attractive Hubbard model, this behaviour is equally expected since 
it is known that the ground state is a Luther-Emery liquid with gapped single-particle 
excitations [8|. Since we observe a continuous crossover from low to high frequencies in 
figure []Jd), we argue that our data supports the existence of a single-particle gap for 
all considered values of Uq] this is confirmed in figure [3] and discussed in the following 
subsection. 

Our analysis of the equal-time correlation functions in this section suggests the 
existence of a crossover from dominant diagonal (charge-density-wave) to dominant off- 
diagonal (pairing) fluctuations as a function of the phonon frequency. This evolution 
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Figure 3. Single-particle spectral function for (a) u)o/t = 0.1 and (b) uio/t = 4. Here 
A = 0.35, n = 0.5, L = 28, and fit = 28. 



is confirmed by the results for P(r) and N(r) = (n r n ) — (^o) shown in figure [21 
where we plot these two correlators at the largest distance r = L/2 for different L. 
Whereas for uo/t = 0.1 the charge correlations decay significantly slower than the pairing 
correlations, see figure EJ^a), dominant pairing correlations can be seen for uo/t = 4 in 
figure[2]^b). The ratio of the exponents obtained from least-square fits is t/cdw/^ss ~ 0.3 
for uo/t = 0.1, and //cdwA/ss ~ 1-4 for u /t = 4. 

3.2. Dynamical correlation functions 

With the help of a stochastic analytical continuation scheme (42], we can extract 
momentum and frequency dependent spectral functions from the imaginary time Green's 
functions accessible in the quantum Monte Carlo simulations. Here we consider the 
single-particle spectral function 

A(k,u) = ^^2(e- pEn +e- pEm ) (11) 

n,m 

x \(n\c ka \m)\ 2 5(E n - E m - u) , 
the dynamical charge structure factor 

N(q,u;) = ^ J £ e ~ PEm \( n \^\ m )\ 2 ( 12 ) 

n,m 

x 6(E n -E m -u), 

where h q = L^^^-e^hj and with the sum rule N(q) = tx^ 1 J duN(q,co), and the 
dynamical spin structure factor S(q, ui) defined as in equation ( fT2l) but with n q replaced 
by S>. 

The results for A(k,oo) shown in figure [3] reveal the existence of a spin gap both at 
low [(a), uo/t = 0.1] and high [(b), uo/t = 4] phonon frequencies, although the size of 
the gap is smaller in figure EJ^b) than in figure [3]Ja). Very similar results have previously 
been obtained for uo/t = 0.1 [T5] . 
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Figure 4. Dynamical spin (a,b) and charge (c,d) structure factors for (a,c) too/t = 0.1 
and (b,d) u /t = 0.5. Here A = 0.35, n = 0.5, L = 28, and (3t = 40. 



The nonzero spin gap is confirmed by our results for the dynamical spin structure 
factor S(q,u) in figure H^a) (for u /t = 0.1) and figure IU(b) (for cu /t = 0.5), which 
reveal the absence of low- lying spectral weight near q = 0. 

The suggested crossover from a Peierls to a superconducting state is expected to 
be reflected in an enhancement of the mobility of the preformed pairs. For the two 
values of the phonon frequency considered here (uo/t = 0.1,0.5), the dynamical charge 
structure factor in figures H](c) and (d) shows two features: a high-energy particle-hole 
continuum reminiscent of the noninteracting problem, and a low-lying excitation band. 
At Uo/t = 0.1, figure 11(c), the dispersion of the latter is extremely flat and within our 
resolution we cannot distinguish the associated charge velocity from zero. This narrow 
mode reflects the very slow dynamics of bipolarons which are predominantly localised 
to form the 2kp charge-density-wave order. At a higher phonon frequency, uj^/t = 0.5, 
see figure S(d), a clear dispersion relation with a minimum at q — 2kp emerges. 

The enhancement of the velocity of the bipolarons as a function of phonon frequency 
is also apparent in the optical conductivity, defined as 

^h^E 6 "^ 1 -^) ( 13 ) 

n,m 

x \(n\j\m)\ 2 5{E n - E m - u) , 
with the paramagnetic current operator j = it ^2i a (c} a c i+1(7 — H.c). Results for o-'(u) 
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Figure 5. Optical conductivity for ujo/t = 0.1 ( ) and tuo/t = 0.5 ( ). Here 

A = 0.35, n = 0.5, L = 28, fit = 40, and a Q = e 2 /2ir. 



are plotted in figure O and reveal that the Drude weight is enhanced by an order of 
magnitude upon increasing the phonon frequency from OJo/t = 0.1 to ooo/t = 0.5. 

The enhancement of the Drude peak follows directly from the dynamical charge 
structure factor. Using the continuity equation, the optical conductivity and the 
dynamical charge structure factor are related via [15] 

a'(q,u) = ^(l-e-^)N( q ,u). (14) 

In the long wavelength limit, we can make the approximation N(q,u) oc q5(u — v c q). 
Here, v c is the charge velocity and the prefactor q stems from phase space available for 
long-wavelength charge flnetnafonsB This ansatz .eads to a\ q = , u) oc v c S(oj), so 
that the enhancement of the velocity v c is reflected in the conductivity 



4. Discussion and conclusions 



We have presented results for static and dynamical correlation functions of the quarter- 
filled Holstein model at A = 0.35, and concentrated on the evolution of the metallic 
state from the adiabatic to the anti-adiabatic limit. For all considered values of the 
phonon frequency, our data support the presence of a spin gap as a consequence of 
backward scattering. Such a system with a gapless charge mode and a gapped spin 
mode is described by the Luther- Emery fixed point. At the latter, the density and 
pairing correlators have the form [8] 

(n(r)n(O)) = ^ + A CO s(2£; F r) (15) 
+ ^k p cos(4/c F r) + . . . , 

( At W A (°)) = ^; + ---' ( 16 ) 

% The factor q equally follows from the sum rule J dujN(q,ui) = N(q) oc q 
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in addition to an exponential decay of the spin-spin correlations. The dominant 
correlations are 2kp density correlations for K p < 1, and pairing correlations for K p > 1. 
Our numerical results, in particular those shown in figure EJ are hence best understood 
in terms of an enhancement of K p from values smaller than one at small oJo/t (leading to 
dominant charge-density-wave correlations) to values larger than 1 in the nonadiabatic 
regime where we observe significantly larger pairing correlations. A quantitative finite- 
size scaling of K p is difficult for several reasons, including the necessity to increase L 
in steps of eight, and the possible importance of logarithmic corrections as a result of 
retardation effects [20] . Nevertheless, a cubic extrapolation based on L = 12,20,28 
gives K p « 0.8 for uj /t = 0.1 and K p « 1.1 for uj /t = 4. 

Since the spin degrees of freedom always remain gapped, we have interpreted 
our results in terms of preformed pairs which are essentially localised in the 2kp 
charge-density- wave state realised in the adiabatic limit. This localisation leads to 
a metallic state with a small Drude weight. Upon increasing the phonon frequency, 
enhanced lattice fluctuations permit bipolaron motion, thereby providing a kinetic 
energy gain. The resulting state is characterised by a significantly larger Drude response 
or, equivalently, an enhanced bipolaron velocity. 

A fruitful direction for future research concerning the evolution of the Peierls state 
is the effect of higher dimensions. By studying, for example, weakly coupled Holstein 
chains, the crossover from a charge-density-wave state to a superfluid discussed here 
may evolve into a phase transition expected to be in the XY universality class. 
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